Univariate Arima Models

Predicting COVID-19 Deaths (STAT 390)

1. Load Libraries, Seed, & Data

Code

# Load necessary libraries
library(forecast)
library(tseries)
library(tidyverse)
library(tidymodels)
library(kableExtra)
library(modeltime)
library(timetk)
library(cowplot)

# setting the seed
set.seed(123)

# Loading Data  ----------------------------------------------------------------
daily_deaths <- read_csv("data/new_daily_deaths.csv") %>% 
  select(daily_deaths)

east <- read_csv("data/east_daily.csv")

midwest <- read_csv("data/midwest_daily.csv")

south <- read_csv("data/south_daily.csv")

west <- read_csv("data/west_daily.csv")

2. Data Splitting

Code

# Data Splitting ---------------------------------------------------------------
east_splits <- time_series_split(east, initial = "1022 days", assess = "114 days")
east_train <- training(east_splits)
east_test <- testing(east_splits)

midwest_splits <- time_series_split(midwest, initial = "1022 days", assess = "114 days")
midwest_train <- training(midwest_splits)
midwest_test <- testing(midwest_splits)

south_splits <- time_series_split(south, initial = "1022 days", assess = "114 days")
south_train <- training(south_splits)
south_test <- testing(south_splits)

west_splits <- time_series_split(west, initial = "1022 days", assess = "114 days")
west_train <- training(west_splits)
west_test <- testing(west_splits)

## Convert Into Time Series ----------------------------------------------------
east_train_ts <- ts(east_train %>%
                      select(daily_deaths))
east_test_ts <- ts(east_test %>%
                     select(daily_deaths))
east_ts <- ts(east %>% 
                select(daily_deaths))

midwest_train_ts <- ts(midwest_train %>%
                      select(daily_deaths))
midwest_test_ts <- ts(midwest_test %>%
                     select(daily_deaths))
midwest_ts <- ts(midwest %>% 
                select(daily_deaths))

south_train_ts <- ts(south_train %>%
                    select(daily_deaths))
south_test_ts <- ts(south_test %>%
                   select(daily_deaths))
south_ts <- ts(south %>%
                 select(daily_deaths))

west_train_ts <- ts(west_train %>%
                   select(daily_deaths))
west_test_ts <- ts(west_test %>%
                  select(daily_deaths))
west_ts <- ts(west %>%
                select(daily_deaths))

## Plotting Training and Testing -----------------------------------------------
par(mfrow=c(2,2))
ggplot() +
  geom_line(data = east_test, aes(x = date, y = daily_deaths, color = "turquoise")) +
  geom_point(data = east_test, aes(x = date, y = daily_deaths, color = "turquoise"), size = .5) +
  geom_line(data = east_train, aes(x = date, y = daily_deaths, color = "orange")) +
  geom_point(data = east_train, aes(x = date, y = daily_deaths, color = "orange"), size = .5) +
  scale_color_manual(name="Data Set", values = c("turquoise", "orange"), labels = c("Training", "Testing")) +
  labs(y = "Daily Deaths",
       x = "Date",
       title = "East Region Training vs. Testing Data")

Code

ggplot() +
  geom_line(data = midwest_test, aes(x = date, y = daily_deaths, color = "turquoise")) +
  geom_point(data = midwest_test, aes(x = date, y = daily_deaths, color = "turquoise"), size = .5) +
  geom_line(data = midwest_train, aes(x = date, y = daily_deaths, color = "orange")) +
  geom_point(data = midwest_train, aes(x = date, y = daily_deaths, color = "orange"), size = .5) +
  scale_color_manual(name="Data Set", values = c("turquoise", "orange"), labels = c("Training", "Testing")) +
  labs(y = "Daily Deaths",
       x = "Date",
       title = "Midwest Region Training vs. Testing Data")

Code

ggplot() +
  geom_line(data = south_test, aes(x = date, y = daily_deaths, color = "turquoise")) +
  geom_point(data = south_test, aes(x = date, y = daily_deaths, color = "turquoise"), size = .5) +
  geom_line(data = south_train, aes(x = date, y = daily_deaths, color = "orange")) +
  geom_point(data = south_train, aes(x = date, y = daily_deaths, color = "orange"), size = .5) +
  scale_color_manual(name="Data Set", values = c("turquoise", "orange"), labels = c("Training", "Testing")) +
  labs(y = "Daily Deaths",
       x = "Date",
       title = "South Region Training vs. Testing Data")

Code

ggplot() +
  geom_line(data = west_test, aes(x = date, y = daily_deaths, color = "turquoise")) +
  geom_point(data = west_test, aes(x = date, y = daily_deaths, color = "turquoise"), size = .5) +
  geom_line(data = west_train, aes(x = date, y = daily_deaths, color = "orange")) +
  geom_point(data = west_train, aes(x = date, y = daily_deaths, color = "orange"), size = .5) +
  scale_color_manual(name="Data Set", values = c("turquoise", "orange"), labels = c("Training", "Testing")) +
  labs(y = "Daily Deaths",
       x = "Date",
       title = "West Region Training vs. Testing Data")

3. Testing for Stationarity

Code

### Using Augmented Dickey-Fuller Test -----------------------------------------
adf_east <- adf.test(east_ts)
print(adf_east) # p-value is 0.01 < 0.05, implying it is stationary

    Augmented Dickey-Fuller Test

data:  east_ts
Dickey-Fuller = -4.5566, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary

Code

adf_midwest <- adf.test(midwest_ts)
print(adf_midwest) # p-value is 0.325 > 0.05, implying it is NOT stationary

    Augmented Dickey-Fuller Test

data:  midwest_ts
Dickey-Fuller = -2.5984, Lag order = 10, p-value = 0.325
alternative hypothesis: stationary

Code

midwest_ts <- diff(midwest_ts)
adf_midwest <- adf.test(midwest_ts)
print(adf_midwest) # p-value is 0.01 < 0.05, implying it is stationary with 1 diff

    Augmented Dickey-Fuller Test

data:  midwest_ts
Dickey-Fuller = -13.43, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary

Code

adf_south <- adf.test(south_ts)
print(adf_south) # p-value is 0.3018 > 0.05, implying it is NOT stationary

    Augmented Dickey-Fuller Test

data:  south_ts
Dickey-Fuller = -2.6533, Lag order = 10, p-value = 0.3018
alternative hypothesis: stationary

Code

south_ts <- diff(south_ts)
adf_south <- adf.test(south_ts)
print(adf_south) # p-value is 0.01 < 0.05, implying it is stationary with 1 diff

    Augmented Dickey-Fuller Test

data:  south_ts
Dickey-Fuller = -11.36, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary

Code

adf_west <- adf.test(west_ts)
print(adf_west) # p-value is 0.4098 > 0.05, implying it is NOT stationary

    Augmented Dickey-Fuller Test

data:  west_ts
Dickey-Fuller = -2.3981, Lag order = 10, p-value = 0.4098
alternative hypothesis: stationary

Code

west_ts <- diff(west_ts)
adf_west <- adf.test(west_ts)
print(adf_west)  # p-value is 0.01 < 0.05, implying it is stationary with 1 diff

    Augmented Dickey-Fuller Test

data:  west_ts
Dickey-Fuller = -12.247, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary

Code

### Making Testing/Training Stationary -----------------------------------------
midwest_train_ts <- diff(midwest_train_ts)
midwest_test_ts <- diff(midwest_test_ts)

midwest_train_ts <- as.data.frame(midwest_train_ts) %>% 
  mutate(daily_deaths = case_when((daily_deaths < 0) ~ 0,
                                  .default = daily_deaths))
midwest_test_ts <- as.data.frame(midwest_test_ts) %>% 
  mutate(daily_deaths = case_when((daily_deaths < 0) ~ 0,
                                  .default = daily_deaths))

midwest_train_ts <- ts(midwest_train_ts)
midwest_test_ts <- ts(midwest_test_ts)

south_train_ts <- diff(south_train_ts)
south_test_ts <- diff(south_test_ts)

south_train_ts <- as.data.frame(south_train_ts) %>% 
  mutate(daily_deaths = case_when((daily_deaths < 0) ~ 0,
                                  .default = daily_deaths))
south_test_ts <- as.data.frame(south_test_ts) %>% 
  mutate(daily_deaths = case_when((daily_deaths < 0) ~ 0,
                                  .default = daily_deaths))

south_train_ts <- ts(south_train_ts)
south_test_ts <- ts(south_test_ts)

west_train_ts <- diff(west_train_ts)
west_test_ts <- diff(west_test_ts)

west_train_ts <- as.data.frame(west_train_ts) %>% 
  mutate(daily_deaths = case_when((daily_deaths < 0) ~ 0,
                                  .default = daily_deaths))
west_test_ts <- as.data.frame(west_test_ts) %>% 
  mutate(daily_deaths = case_when((daily_deaths < 0) ~ 0,
                                  .default = daily_deaths))

west_train_ts <- ts(west_train_ts)
west_test_ts <- ts(west_test_ts)

4. ACF & PACF Plots

Code

## ACF & PACF Plots ------------------------------------------------------------
par(mfrow=c(3, 1))
acf(east_ts, main = "East")
pacf(east_ts, main = "East")
plot(east_ts, main = "Daily Deaths: East")

Code

acf(midwest_ts, main = "Midwest")
pacf(midwest_ts, main = "Midwest")
plot(midwest_ts, main = "Daily Deaths: Midwest")

Code

acf(south_ts, main = "South")
pacf(south_ts, main = "South")
plot(south_ts, main = "Daily Deaths: South")

Code

acf(west_ts, main = "West")
pacf(west_ts, main = "West")
plot(south_ts, main = "Daily Deaths: South")

Code

plot(west_ts, main = "Daily Deaths: West")

5. Time Series Decomposition

Code

# total daily death decomposition
daily_ts <- ts(daily_deaths, frequency = 365)

decomposed_ts <- decompose(daily_ts)

# Plot the original time series and the decomposed components
par(mfrow=c(2,1))
plot(daily_ts, main = "Total Original Time Series", ylab = "Value", col = "blue")
plot(decomposed_ts$trend, main = "Trend Component", ylab = "Value", col = "red")

Code

plot(decomposed_ts$seasonal, main = "Seasonal Component", ylab = "Value", col = "green")
plot(decomposed_ts$random, main = "Residual Component", ylab = "Value", col = "purple")

Code

# East Decomposition -----------------------------------------------------------
east_daily <- east %>% 
  select(daily_deaths)

east_ts <- ts(east_daily, frequency = 365)

decomposed_east <- decompose(east_ts)

# Plot the original time series and the decomposed components
par(mfrow=c(2,1))
plot(east_ts, main = "East Original Time Series", ylab = "Value", col = "blue")
plot(decomposed_east$trend, main = "Trend Component", ylab = "Value", col = "red")

Code

plot(decomposed_east$seasonal, main = "Seasonal Component", ylab = "Value", col = "green")
plot(decomposed_east$random, main = "Residual Component", ylab = "Value", col = "purple")

Code

## Midwest Decomposition ----------------------------------------------------------
midwest_daily <- midwest %>% 
  select(daily_deaths)

midwest_ts <- ts(midwest_daily, frequency = 365)

decomposed_midwest <- decompose(midwest_ts)

# Plot the original time series and the decomposed components
par(mfrow=c(2,1))
plot(midwest_ts, main = "Midwest Original Time Series", ylab = "Value", col = "blue")
plot(decomposed_midwest$trend, main = "Trend Component", ylab = "Value", col = "red")

Code

plot(decomposed_midwest$seasonal, main = "Seasonal Component", ylab = "Value", col = "green")
plot(decomposed_midwest$random, main = "Residual Component", ylab = "Value", col = "purple")

Code

## South Decomposition ----------------------------------------------------------
south_daily <- south %>% 
  select(daily_deaths)

south_ts <- ts(south_daily, frequency = 365)

decomposed_south <- decompose(south_ts)

# Plot the original time series and the decomposed components
par(mfrow=c(2,1))
plot(south_ts, main = "South Original Time Series", ylab = "Value", col = "blue")
plot(decomposed_south$trend, main = "Trend Component", ylab = "Value", col = "red")

Code

plot(decomposed_south$seasonal, main = "Seasonal Component", ylab = "Value", col = "green")
plot(decomposed_south$random, main = "Residual Component", ylab = "Value", col = "purple")

Code

## West Decomposition ----------------------------------------------------------
west_daily <- west %>% 
  select(daily_deaths)

west_ts <- ts(west_daily, frequency = 365)

decomposed_west <- decompose(west_ts)

# Plot the original time series and the decomposed components
par(mfrow=c(2,1))
plot(west_ts, main = "West Original Time Series", ylab = "Value", col = "blue")
plot(decomposed_west$trend, main = "Trend Component", ylab = "Value", col = "red")

Code

plot(decomposed_west$seasonal, main = "Seasonal Component", ylab = "Value", col = "green")
plot(decomposed_west$random, main = "Residual Component", ylab = "Value", col = "purple")

6. ARIMA Model

6.1 Initial Model

Code

# Building ARIMA Model ---------------------------------------------------------
east_arima <- arima(east_train_ts, order=c(1,0,1))
midwest_arima <- arima(midwest_train_ts, order=c(1,1,1))
south_arima <- arima(south_train_ts, order=c(1,1,1))
west_arima <- arima(west_train_ts, order=c(1,1,1))

## Summary & Forecast ---------------------------------------------------------
summary(east_arima)

Call:
arima(x = east_train_ts, order = c(1, 0, 1))

Coefficients:
         ar1      ma1  intercept
      0.9176  -0.0388   197.3738
s.e.  0.0147   0.0493    39.7998

sigma^2 estimated as 12153:  log likelihood = -6257.17,  aic = 12522.33

Training set error measures:
                    ME     RMSE      MAE  MPE MAPE      MASE        ACF1
Training set 0.1151898 110.2407 64.66739 -Inf  Inf 0.9918341 0.008515868

Code

summary(midwest_arima)

Call:
arima(x = midwest_train_ts, order = c(1, 1, 1))

Coefficients:
          ar1      ma1
      -0.2237  -0.9343
s.e.   0.0313   0.0090

sigma^2 estimated as 16650:  log likelihood = -6405.86,  aic = 12817.73

Training set error measures:
                    ME     RMSE     MAE MPE MAPE      MASE        ACF1
Training set 0.7797006 128.9735 78.3255 NaN  Inf 0.7410171 -0.06964661

Code

summary(south_arima)

Call:
arima(x = south_train_ts, order = c(1, 1, 1))

Coefficients:
          ar1      ma1
      -0.1913  -0.9453
s.e.   0.0317   0.0095

sigma^2 estimated as 44608:  log likelihood = -6908.51,  aic = 13823.02

Training set error measures:
                   ME     RMSE     MAE MPE MAPE      MASE        ACF1
Training set 1.277928 211.1023 129.556 NaN  Inf 0.7498262 -0.01884038

Code

summary(west_arima)

Call:
arima(x = west_train_ts, order = c(1, 1, 1))

Coefficients:
          ar1      ma1
      -0.0842  -0.9603
s.e.   0.0319   0.0072

sigma^2 estimated as 10324:  log likelihood = -6162.22,  aic = 12330.44

Training set error measures:
                    ME    RMSE      MAE  MPE MAPE      MASE         ACF1
Training set 0.6807731 101.558 62.19197 -Inf  Inf 0.7886397 -0.009430079

Code

east_values <- forecast(east_arima, h = length(east_test_ts))
midwest_values <- forecast(midwest_arima, h = length(midwest_test_ts))
south_values <- forecast(south_arima, h = length(south_test_ts))
west_values <- forecast(west_arima, h = length(west_test_ts))

6.2 Accuracy

Code

## Evaluating Accuracy ---------------------------------------------------------
east_initial_arima_acc <- forecast::accuracy(east_values)

east_initial_arima_acc %>% 
  kbl(caption = "<center>Initial East Accuracy<center>") %>% 
  kable_styling()

Initial East Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.1151898	110.2407	64.66739	-Inf	Inf	0.9918341	0.0085159

Code

### East Region
east_test_plot <- as.data.frame(east_test_ts) %>% 
  mutate(num = seq(1023, (1023+113)))

east_initial_arima_plot <- forecast::autoplot(east_values) +
  geom_line(data = east_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "East Region Initial ARIMA",
       y = "Daily Deaths")

east_initial_arima_plot

Code

### Midwest Region

midwest_initial_arima_acc <- forecast::accuracy(midwest_values)

midwest_initial_arima_acc %>% 
  kbl(caption = "<center>Initial Midwest Accuracy<center>") %>% 
  kable_styling()

Initial Midwest Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.7797006	128.9735	78.3255	NaN	Inf	0.7410171	-0.0696466

Code

midwest_test_plot <- as.data.frame(midwest_test_ts) %>% 
  mutate(num = seq(1023, (1023+112)))

midwest_initial_arima_plot <- forecast::autoplot(midwest_values) +
  geom_line(data = midwest_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "Midwest Region Initial ARIMA",
       y = "Daily Deaths")

midwest_initial_arima_plot

Code

### South Region

south_initial_arima_acc <- forecast::accuracy(south_values)

south_initial_arima_acc %>% 
  kbl(caption = "<center>Initial South Accuracy<center>") %>% 
  kable_styling()

Initial South Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	1.277928	211.1023	129.556	NaN	Inf	0.7498262	-0.0188404

Code

south_test_plot <- as.data.frame(south_test_ts) %>% 
  mutate(num = seq(1023, (1023+112)))

south_initial_arima_plot <- forecast::autoplot(south_values) +
  geom_line(data = south_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "South Region Initial ARIMA",
       y = "Daily Deaths")

south_initial_arima_plot 

Code

### West Region

west_initial_arima_acc <- forecast::accuracy(west_values)

west_initial_arima_acc %>% 
  kbl(caption = "<center>Initial West Accuracy<center>") %>% 
  kable_styling()

Initial West Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.6807731	101.558	62.19197	-Inf	Inf	0.7886397	-0.0094301

Code

west_test_plot <- as.data.frame(west_test_ts) %>% 
  mutate(num = seq(1023, (1023+112)))

west_initial_arima_plot <- forecast::autoplot(west_values) +
  geom_line(data = west_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "West Region Initial ARIMA",
       y = "Daily Deaths") 

west_initial_arima_plot

6.3 Grid Search

Code

# Grid Search for p and q ------------------------------------------------------

p_values <- 0:5  # AR order
d_values <- 0:0  # I order
q_values <- 0:5  # MA order

## East - Perform grid search for ARIMA parameters  ----------------------------
best_east_model <- NULL
best_east_aic <- Inf

for (p in p_values) {
  for (d in d_values) {
    for (q in q_values) {
      # Fit ARIMA model
      arima_model <- tryCatch(
        {
          fit <- Arima(east_ts, order = c(p, d, q))
          fit
        },
        error = function(e) {
          NULL
        }
      )
      
      # Check if ARIMA model was successfully fitted
      if (!is.null(arima_model)) {
        # Check AIC criterion
        if (AIC(arima_model) < best_east_aic) {
          best_east_model <- arima_model
          best_east_aic <- AIC(arima_model)
          best_east_order <- c(p, d, q)
        }
      }
    }
  }
}

Code

## Midwest - Perform grid search for ARIMA parameters  -------------------------
p_values <- 0:5  # AR order
d_values <- 0:1  # I order
q_values <- 0:5  # MA order

best_midwest_model <- NULL
best_midwest_aic <- Inf

for (p in p_values) {
  for (d in d_values) {
    for (q in q_values) {
      # Fit ARIMA model
      arima_model <- tryCatch(
        {
          fit <- Arima(midwest_ts, order = c(p, d, q))
          fit
        },
        error = function(e) {
          NULL
        }
      )
      
      # Check if ARIMA model was successfully fitted
      if (!is.null(arima_model)) {
        # Check AIC criterion
        if (AIC(arima_model) < best_midwest_aic) {
          best_midwest_model <- arima_model
          best_midwest_aic <- AIC(arima_model)
          best_midwest_order <- c(p, d, q)
        }
      }
    }
  }
}

Code

## South - Perform grid search for ARIMA parameters  ---------------------------
p_values <- 0:5  # AR order
d_values <- 0:1  # I order
q_values <- 0:5  # MA order

best_south_model <- NULL
best_south_aic <- Inf

for (p in p_values) {
  for (d in d_values) {
    for (q in q_values) {
      # Fit ARIMA model
      arima_model <- tryCatch(
        {
          fit <- Arima(south_ts, order = c(p, d, q))
          fit
        },
        error = function(e) {
          NULL
        }
      )
      
      # Check if ARIMA model was successfully fitted
      if (!is.null(arima_model)) {
        # Check AIC criterion
        if (AIC(arima_model) < best_south_aic) {
          best_south_model <- arima_model
          best_south_aic <- AIC(arima_model)
          best_south_order <- c(p, d, q)
        }
      }
    }
  }
}

Code

## West - Perform grid search for ARIMA parameters  ----------------------------
p_values <- 0:5  # AR order
d_values <- 0:1  # I order
q_values <- 0:5  # MA order

best_west_model <- NULL
best_west_aic <- Inf

for (p in p_values) {
  for (d in d_values) {
    for (q in q_values) {
      # Fit ARIMA model
      arima_model <- tryCatch(
        {
          fit <- Arima(west_ts, order = c(p, d, q))
          fit
        },
        error = function(e) {
          NULL
        }
      )
      
      # Check if ARIMA model was successfully fitted
      if (!is.null(arima_model)) {
        # Check AIC criterion
        if (AIC(arima_model) < best_west_aic) {
          best_west_model <- arima_model
          best_west_aic <- AIC(arima_model)
          best_west_order <- c(p, d, q)
        }
      }
    }
  }
}

Code

# Saving results
save(best_east_order, best_east_model,
     best_south_order, best_south_model,
     best_midwest_order, best_midwest_model,
     best_west_order, best_west_model, 
     file = "results/arima_grid.rda")

Best East Model

Code

# load results
load("results/arima_grid.rda")
# Display the best ARIMA model and its parameters
print(paste(best_east_order))

[1] "5" "0" "5"

Code

print(best_east_model)

Series: east_ts 
ARIMA(5,0,5) with non-zero mean 

Coefficients:
         ar1      ar2     ar3     ar4      ar5      ma1     ma2     ma3
      2.5325  -2.4820  0.7383  0.6063  -0.4009  -2.1286  1.6389  0.0590
s.e.  0.1056   0.3215  0.4321  0.3068   0.0949   0.0975  0.2582  0.3121
          ma4     ma5      mean
      -0.8826  0.4627  184.1912
s.e.   0.1987  0.0577   60.5578

sigma^2 = 6735:  log likelihood = -6616.45
AIC=13256.9   AICc=13257.17   BIC=13317.32

Best Midwest Model

Code

print(paste(best_midwest_order))

[1] "5" "1" "5"

Code

print(best_midwest_model)

Series: midwest_ts 
ARIMA(5,1,5) 

Coefficients:
         ar1      ar2      ar3      ar4      ar5      ma1     ma2      ma3
      0.2188  -0.9672  -0.0252  -0.5260  -0.5481  -1.1980  1.4606  -1.2997
s.e.  0.0353   0.0286   0.0485   0.0288   0.0324   0.0399  0.0491   0.0522
         ma4      ma5
      1.0143  -0.4534
s.e.  0.0483   0.0337

sigma^2 = 15644:  log likelihood = -7089.78
AIC=14201.56   AICc=14201.8   BIC=14256.94

Best South Model

Code

print(paste(best_south_order))

[1] "5" "1" "5"

Code

print(best_south_model)

Series: south_ts 
ARIMA(5,1,5) 

Coefficients:
         ar1      ar2     ar3      ar4      ar5      ma1     ma2      ma3
      0.4802  -1.1648  0.3304  -0.7153  -0.3031  -1.3248  1.5824  -1.3808
s.e.  0.0544   0.0456  0.0758   0.0423   0.0488   0.0519  0.0749   0.0844
         ma4      ma5
      1.0487  -0.4103
s.e.  0.0652   0.0413

sigma^2 = 38413:  log likelihood = -7598.92
AIC=15219.84   AICc=15220.08   BIC=15275.22

Best West Model

Code

print(paste(best_west_order))

[1] "5" "1" "5"

Code

print(best_west_model)

Series: west_ts 
ARIMA(5,1,5) 

Coefficients:
         ar1      ar2     ar3      ar4      ar5      ma1     ma2      ma3
      0.4629  -1.1476  0.3025  -0.7075  -0.3255  -1.3275  1.6081  -1.4132
s.e.  0.0516   0.0452  0.0734   0.0403   0.0465   0.0492  0.0746   0.0849
         ma4      ma5
      1.1110  -0.4753
s.e.  0.0595   0.0368

sigma^2 = 8126:  log likelihood = -6717.71
AIC=13457.43   AICc=13457.66   BIC=13512.8

6.4 Rebuilding ARIMA Model

Code

# Rebuilding ARIMA Models ------------------------------------------------------
east_arima <- arima(east_train_ts, order=c(5,0,5))
midwest_arima <- arima(midwest_train_ts, order=c(5,1,5))
south_arima <- arima(south_train_ts, order=c(5,1,5))
west_arima <- arima(west_train_ts, order=c(5,1,5), method = "ML")


east_values <- forecast(east_arima, h = length(east_test_ts))
midwest_values <- forecast(midwest_arima, h = length(midwest_test_ts))
south_values <- forecast(south_arima, h = length(south_test_ts))
west_values <- forecast(west_arima, h = length(west_test_ts))

## Evaluating Accuracy & Plots -------------------------------------------------

east_tuned_arima_acc <- forecast::accuracy(east_values)

east_tuned_arima_acc %>% 
  kbl(caption = "<center>Tuned East Accuracy<center>") %>% 
  kable_styling()

Tuned East Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.1095068	80.35755	47.18823	NaN	Inf	0.723748	-0.0254917

Code

### East Region
east_tuned_arima_plot <- forecast::autoplot(east_values) +
  geom_line(data = east_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "East Region Tuned ARIMA",
       y = "Daily Deaths")

east_tuned_arima_plot

Code

### Midwest Region

midwest_tuned_arima_acc <- forecast::accuracy(midwest_values)

midwest_tuned_arima_acc %>% 
  kbl(caption = "<center>Tuned Midwest Accuracy<center>") %>% 
  kable_styling()

Tuned Midwest Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.4882952	96.70073	54.43655	NaN	Inf	0.5150099	-0.0118025

Code

midwest_tuned_arima_plot <- forecast::autoplot(midwest_values) +
  geom_line(data = midwest_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "Midwest Region Tuned ARIMA",
       y = "Daily Deaths")

midwest_tuned_arima_plot

Code

### South Region

south_tuned_arima_acc <- forecast::accuracy(south_values)

south_tuned_arima_acc %>% 
  kbl(caption = "<center>Tuned South Accuracy<center>") %>% 
  kable_styling()

Tuned South Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.7086084	172.0968	94.70831	NaN	Inf	0.5481396	-0.0036251

Code

south_tuned_arima_plot <- forecast::autoplot(south_values) +
  geom_line(data = south_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "South Region Tuned ARIMA",
       y = "Daily Deaths")

south_tuned_arima_plot 

Code

### West Region
west_tuned_arima_acc <- forecast::accuracy(west_values)

west_tuned_arima_acc %>% 
  kbl(caption = "<center>Tuned West Accuracy<center>") %>% 
  kable_styling()

Tuned West Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.4533154	79.29443	46.88033	NaN	Inf	0.5944769	-0.0091414

Code

west_tuned_arima_plot <- forecast::autoplot(west_values) +
  geom_line(data = west_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "West Region Tuned ARIMA",
       y = "Daily Deaths") 

west_tuned_arima_plot

Code

# saving results
save(east_tuned_arima_acc, east_initial_arima_acc, east_initial_arima_plot, east_tuned_arima_plot,
     midwest_tuned_arima_acc, midwest_initial_arima_acc, midwest_initial_arima_plot, midwest_tuned_arima_plot,
     south_tuned_arima_acc, south_initial_arima_acc, south_initial_arima_plot, south_tuned_arima_plot,
     west_tuned_arima_acc, west_initial_arima_acc, west_initial_arima_plot, west_tuned_arima_plot,
     file = "results/uni_arima.rda"
     )

7. SARIMA Model

7.1 Initial SARIMA Model

Code

# Building SARIMA model  ------------------------------------------------------
east_sarima <- arima(east_train_ts, seasonal = list(order = c(1,1,1)))
midwest_sarima <- arima(midwest_train_ts, seasonal = list(order = c(1,1,1)))
south_sarima <- arima(south_train_ts, seasonal = list(order = c(1,1,1)))
west_sarima <- arima(west_train_ts, seasonal = list(order = c(1,1,1)))

# Summary of the SARIMA model
summary(east_sarima)

Call:
arima(x = east_train_ts, seasonal = list(order = c(1, 1, 1)))

Coefficients:
        sar1     sma1
      0.4706  -0.8016
s.e.  0.0382   0.0211

sigma^2 estimated as 11156:  log likelihood = -6206.64,  aic = 12419.27

Training set error measures:
                    ME     RMSE     MAE  MPE MAPE      MASE      ACF1
Training set 0.1606608 105.5715 61.2767 -Inf  Inf 0.9398295 0.1140208

Code

summary(midwest_sarima)

Call:
arima(x = midwest_train_ts, seasonal = list(order = c(1, 1, 1)))

Coefficients:
         sar1     sma1
      -0.2237  -0.9343
s.e.   0.0313   0.0090

sigma^2 estimated as 16650:  log likelihood = -6405.86,  aic = 12817.73

Training set error measures:
                    ME     RMSE     MAE MPE MAPE      MASE        ACF1
Training set 0.7797006 128.9735 78.3255 NaN  Inf 0.7410171 -0.06964661

Code

summary(south_sarima)

Call:
arima(x = south_train_ts, seasonal = list(order = c(1, 1, 1)))

Coefficients:
         sar1     sma1
      -0.1913  -0.9453
s.e.   0.0317   0.0095

sigma^2 estimated as 44608:  log likelihood = -6908.51,  aic = 13823.02

Training set error measures:
                   ME     RMSE     MAE MPE MAPE      MASE        ACF1
Training set 1.277928 211.1023 129.556 NaN  Inf 0.7498262 -0.01884038

Code

summary(west_sarima)

Call:
arima(x = west_train_ts, seasonal = list(order = c(1, 1, 1)))

Coefficients:
         sar1     sma1
      -0.0842  -0.9603
s.e.   0.0319   0.0072

sigma^2 estimated as 10324:  log likelihood = -6162.22,  aic = 12330.44

Training set error measures:
                    ME    RMSE      MAE  MPE MAPE      MASE         ACF1
Training set 0.6807731 101.558 62.19197 -Inf  Inf 0.7886397 -0.009430079

Code

## Forecasting with the SARIMA model -------------------------------------------
east_values <- forecast(east_sarima, h = length(east_test_ts))
midwest_values <- forecast(midwest_sarima, h = length(midwest_test_ts))
south_values <- forecast(south_sarima, h = length(south_test_ts))
west_values <- forecast(west_sarima, h = length(west_test_ts))

7.2 Accuracy

Code

## Plotting & Accuracy

## East Region
east_initial_sarima_acc <- forecast::accuracy(east_values)

east_initial_sarima_acc %>% 
  kbl(caption = "<center>Initial SARIMA East Accuracy<center>") %>% 
  kable_styling()

Initial SARIMA East Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.1606608	105.5715	61.2767	-Inf	Inf	0.9398295	0.1140208

Code

east_initial_sarima_plot <- forecast::autoplot(east_values) +
  geom_line(data = east_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "East Region Initial SARIMA",
       y = "Daily Deaths")

east_initial_sarima_plot

Code

## Midwest Region

midwest_initial_sarima_acc <- forecast::accuracy(midwest_values)

midwest_initial_sarima_acc %>% 
  kbl(caption = "<center>Initial SARIMA Midwest Accuracy<center>") %>% 
  kable_styling()

Initial SARIMA Midwest Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.7797006	128.9735	78.3255	NaN	Inf	0.7410171	-0.0696466

Code

midwest_initial_sarima_plot <- forecast::autoplot(midwest_values) +
  geom_line(data = midwest_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "Midwest Region Initial SARIMA",
       y = "Daily Deaths")

midwest_initial_sarima_plot

Code

## South Region
south_initial_sarima_acc <- forecast::accuracy(south_values)

south_initial_sarima_acc %>%
  kbl(caption = "<center>Initial SARIMA South Accuracy<center>") %>%
  kable_styling()

Initial SARIMA South Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	1.277928	211.1023	129.556	NaN	Inf	0.7498262	-0.0188404

Code

south_initial_sarima_plot <- forecast::autoplot(south_values) +
  geom_line(data = south_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "South Region Initial SARIMA",
       y = "Daily Deaths")

south_initial_sarima_plot

Code

## West Region

west_initial_sarima_acc <- forecast::accuracy(west_values)

west_initial_sarima_acc %>% 
  kbl(caption = "<center>Initial SARIMA West Accuracy<center>") %>% 
  kable_styling()

Initial SARIMA West Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.6807731	101.558	62.19197	-Inf	Inf	0.7886397	-0.0094301

Code

west_initial_sarima_plot <- forecast::autoplot(west_values) +
  geom_line(data = west_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "West Region Initial SARIMA",
       y = "Daily Deaths")

west_initial_sarima_plot

7.3 Grid Search

Code

## Define the parameter grids --------------------------------------------------
p_grid <- 0:5  # AR parameter
d_grid <- 0:0  # Differencing parameter
q_grid <- 0:5  # MA parameter
P_grid <- 0:1  # Seasonal AR parameter
D_grid <- 0:1  # Seasonal differencing parameter
Q_grid <- 0:1  # Seasonal MA parameter
s <- 12        # Seasonal period



## East Grid Search ------------------------------------------------------------

# Create East Region
east_train <- east_train %>% 
  select(daily_deaths)

# Create an empty data frame to store the results
east_results <- data.frame(order = character(),
                      seasonal = character(),
                      AIC = numeric(),
                      stringsAsFactors = FALSE)

# Perform grid search
for (p in p_grid) {
  for (d in d_grid) {
    for (q in q_grid) {
      for (P in P_grid) {
        for (D in D_grid) {
          for (Q in Q_grid) {
            order <- c(p, d, q)
            seasonal <- c(P, D, Q, s)
            model <- tryCatch(
              {
                fit <- arima(east_train, order = order, seasonal = seasonal)
                AIC_val <- AIC(fit)
                # Append the results to the data frame
                east_results <- rbind(east_results, data.frame(order = paste(order, collapse = ", "),
                                                     seasonal = paste(seasonal, collapse = ", "),
                                                     AIC = AIC_val))
              },
              error = function(e) NULL
            )
          }
        }
      }
    }
  }
}

d_grid <- 0:1  # Differencing parameter

## Midwest Grid Search ------------------------------------------------------------

# Create Midwest Region
midwest_train <- midwest_train %>% 
  select(daily_deaths)

# Create an empty data frame to store the results
midwest_results <- data.frame(order = character(),
                           seasonal = character(),
                           AIC = numeric(),
                           stringsAsFactors = FALSE)

# Perform grid search
for (p in p_grid) {
  for (d in d_grid) {
    for (q in q_grid) {
      for (P in P_grid) {
        for (D in D_grid) {
          for (Q in Q_grid) {
            order <- c(p, d, q)
            seasonal <- c(P, D, Q, s)
            model <- tryCatch(
              {
                fit <- arima(midwest_train, order = order, seasonal = seasonal)
                AIC_val <- AIC(fit)
                # Append the results to the data frame
                midwest_results <- rbind(midwest_results, data.frame(order = paste(order, collapse = ", "),
                                                               seasonal = paste(seasonal, collapse = ", "),
                                                               AIC = AIC_val))
              },
              error = function(e) NULL
            )
          }
        }
      }
    }
  }
}

## South Grid Search ------------------------------------------------------------

# Create South Region
south_train <- south_train %>% 
  select(daily_deaths)

# Create an empty data frame to store the results
south_results <- data.frame(order = character(),
                           seasonal = character(),
                           AIC = numeric(),
                           stringsAsFactors = FALSE)

# Perform grid search
for (p in p_grid) {
  for (d in d_grid) {
    for (q in q_grid) {
      for (P in P_grid) {
        for (D in D_grid) {
          for (Q in Q_grid) {
            order <- c(p, d, q)
            seasonal <- c(P, D, Q, s)
            model <- tryCatch(
              {
                fit <- arima(south_train, order = order, seasonal = seasonal)
                AIC_val <- AIC(fit)
                # Append the results to the data frame
                south_results <- rbind(south_results, data.frame(order = paste(order, collapse = ", "),
                                                               seasonal = paste(seasonal, collapse = ", "),
                                                               AIC = AIC_val))
              },
              error = function(e) NULL
            )
          }
        }
      }
    }
  }
}


## West Grid Search ------------------------------------------------------------

# Create West Region
west_train <- west_train %>% 
  select(daily_deaths)

# Create an empty data frame to store the results
west_results <- data.frame(order = character(),
                           seasonal = character(),
                           AIC = numeric(),
                           stringsAsFactors = FALSE)

# Perform grid search
for (p in p_grid) {
  for (d in d_grid) {
    for (q in q_grid) {
      for (P in P_grid) {
        for (D in D_grid) {
          for (Q in Q_grid) {
            order <- c(p, d, q)
            seasonal <- c(P, D, Q, s)
            model <- tryCatch(
              {
                fit <- arima(west_train, order = order, seasonal = seasonal)
                AIC_val <- AIC(fit)
                # Append the results to the data frame
                west_results <- rbind(west_results, data.frame(order = paste(order, collapse = ", "),
                                                               seasonal = paste(seasonal, collapse = ", "),
                                                               AIC = AIC_val))
              },
              error = function(e) NULL
            )
          }
        }
      }
    }
  }
}

Code

save(east_results, midwest_results, south_results, west_results, file = "results/sarima_grid.rda")

Code

load("results/sarima_grid.rda")

Best East SARIMA

Code

# Find the best model (minimum AIC)
best_east_model <- east_results[which.min(east_results$AIC), ]
print(best_east_model)

      order    seasonal      AIC
213 5, 0, 4 0, 1, 1, 12 11865.53

Best Midwest SARIMA

Code

# Midwest Results
# Find the best model (minimum AIC)
best_midwest_model <- midwest_results[which.min(midwest_results$AIC), ]
print(best_midwest_model)

      order    seasonal      AIC
385 5, 1, 5 0, 0, 1, 12 12819.82

Best South SARIMA

Code

# South Results
# Find the best model (minimum AIC)
best_south_model <- south_results[which.min(south_results$AIC), ]
print(best_south_model)

      order    seasonal      AIC
402 5, 1, 5 0, 0, 0, 12 13698.93

Best West SARIMA

Code

# Find the best model (minimum AIC)
best_west_model <- west_results[which.min(west_results$AIC), ]
print(best_west_model)

      order    seasonal      AIC
385 5, 0, 5 1, 1, 0, 12 12141.63

7.4 Rebuilding SARIMA

Code

# Rebuilding SARIMA Models -----------------------------------------------------
east_sarima <- arima(east_train_ts, order = c(5, 0, 4), 
                     seasonal = list(order = c(0, 1, 1), period = 12))
midwest_sarima <- arima(midwest_train_ts, order = c(5, 1, 5),
                        seasonal = list(order = c(0, 0, 1), period = 12))
south_sarima <- arima(south_train_ts, order = c(5, 1, 5),
                           seasonal = list(order = c(0, 0, 0), period = 12))
west_sarima <- arima(west_train_ts, order = c(5, 0, 5),
                     seasonal = list(order = c(1, 1, 0), period = 12))


## Forecasting with the NEW SARIMA model ---------------------------------------
east_values <- forecast(east_sarima, h = length(east_test_ts))
midwest_values <- forecast(midwest_sarima, h = length(midwest_test_ts))
south_values <- forecast(south_sarima, h = length(south_test_ts))
west_values <- forecast(west_sarima, h = length(west_test_ts))

Code

## Plotting & Accuracy

## East Region
east_tuned_sarima_acc <- forecast::accuracy(east_values)

east_tuned_sarima_acc %>% 
  kbl(caption = "<center>Tuned SARIMA East Accuracy<center>") %>% 
  kable_styling()

Tuned SARIMA East Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.1031443	87.81972	53.94287	NaN	Inf	0.8273471	-0.0998164

Code

east_tuned_sarima_plot <- forecast::autoplot(east_values) +
  geom_line(data = east_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "East Region Tuned SARIMA",
       y = "Daily Deaths")

east_tuned_sarima_plot

Code

## Midwest Region
midwest_tuned_sarima_acc <- forecast::accuracy(midwest_values)

midwest_tuned_sarima_acc %>% 
  kbl(caption = "<center>Tuned SARIMA Midwest Accuracy<center>") %>% 
  kable_styling()

Tuned SARIMA Midwest Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.4955934	106.0021	63.95655	NaN	Inf	0.6050762	0.0305554

Code

midwest_tuned_sarima_plot <- forecast::autoplot(midwest_values) +
  geom_line(data = midwest_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "Midwest Region Tuned SARIMA",
       y = "Daily Deaths")

midwest_tuned_sarima_plot

Code

## South Region
south_tuned_sarima_acc <- forecast::accuracy(south_values)

south_tuned_sarima_acc %>% 
  kbl(caption = "<center>Tuned SARIMA South Accuracy<center>") %>% 
  kable_styling()

Tuned SARIMA South Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.7086084	172.0968	94.70831	NaN	Inf	0.5481396	-0.0036251

Code

south_tuned_sarima_plot <- forecast::autoplot(south_values) +
  geom_line(data = south_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "South Region Tuned SARIMA",
       y = "Daily Deaths")

south_tuned_sarima_plot

Code

## West Region
west_tuned_sarima_acc <- forecast::accuracy(west_values)

west_tuned_sarima_acc %>% 
  kbl(caption = "<center>Tuned SARIMA West Accuracy<center>") %>% 
  kable_styling()

Tuned SARIMA West Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.681019	89.05196	55.78604	NaN	Inf	0.7074078	0.0029928

Code

west_tuned_sarima_plot <- forecast::autoplot(west_values) +
  geom_line(data = west_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "West Region Tuned SARIMA",
       y = "Daily Deaths")

west_tuned_sarima_plot

Code

# saving results
save(east_tuned_sarima_acc, east_initial_sarima_acc, east_initial_sarima_plot, east_tuned_sarima_plot,
     midwest_tuned_sarima_acc, midwest_initial_sarima_acc, midwest_initial_sarima_plot, midwest_tuned_sarima_plot,
     south_tuned_sarima_acc, south_initial_sarima_acc, south_initial_sarima_plot, south_tuned_sarima_plot,
     west_tuned_sarima_acc, west_initial_sarima_acc, west_initial_sarima_plot, west_tuned_sarima_plot,
     file = "results/uni_sarima.rda"
     )

8. Auto ARIMA Model

8.1 Building AutoARIMA

Code

# Building AutoARIMA Model -----------------------------------------------------
east_arima <- auto.arima(east_train_ts, seasonal = TRUE)
midwest_arima <- auto.arima(midwest_train_ts, seasonal = TRUE)
south_arima <- auto.arima(south_train_ts, seasonal = TRUE)
west_arima <- auto.arima(west_train_ts, seasonal = TRUE)

# Summary of the ARIMA model
summary(east_arima)

Series: east_train_ts 
ARIMA(5,1,2) 

Coefficients:
         ar1      ar2      ar3      ar4      ar5      ma1     ma2
      0.5142  -0.7304  -0.1394  -0.3479  -0.2207  -1.0295  0.7353
s.e.  0.0480   0.0344   0.0426   0.0329   0.0418   0.0410  0.0265

sigma^2 = 6861:  log likelihood = -5956.39
AIC=11928.79   AICc=11928.93   BIC=11968.22

Training set error measures:
                    ME     RMSE     MAE MPE MAPE      MASE       ACF1
Training set 0.1506172 82.50546 47.1043 NaN  Inf 0.7224607 0.01707355

Code

summary(midwest_arima)

Series: midwest_train_ts 
ARIMA(5,1,4) 

Coefficients:
          ar1      ar2      ar3      ar4      ar5      ma1     ma2      ma3
      -0.0349  -0.7605  -0.0404  -0.1320  -0.4194  -1.1696  0.9929  -0.9722
s.e.   0.0461   0.0370   0.0572   0.0371   0.0377   0.0441  0.0492   0.0459
         ma4
      0.4244
s.e.  0.0381

sigma^2 = 10525:  log likelihood = -6168.66
AIC=12357.32   AICc=12357.54   BIC=12406.6

Training set error measures:
                   ME     RMSE      MAE MPE MAPE      MASE        ACF1
Training set 0.385353 102.0896 59.56282 NaN  Inf 0.5635082 -0.05532994

Code

summary(south_arima)

Series: south_train_ts 
ARIMA(1,1,2) 

Coefficients:
         ar1      ma1     ma2
      0.4799  -1.6879  0.7236
s.e.  0.0578   0.0428  0.0405

sigma^2 = 42636:  log likelihood = -6884.08
AIC=13776.16   AICc=13776.2   BIC=13795.87

Training set error measures:
                    ME     RMSE      MAE MPE MAPE      MASE        ACF1
Training set 0.7155848 206.0793 126.7407 NaN  Inf 0.7335323 0.008694211

Code

summary(west_arima)

Series: west_train_ts 
ARIMA(5,1,3) 

Coefficients:
         ar1      ar2      ar3      ar4      ar5      ma1     ma2      ma3
      0.4572  -0.6143  -0.2778  -0.2654  -0.2584  -1.8532  1.6316  -0.6651
s.e.  0.0588   0.0355   0.0387   0.0340   0.0435   0.0544  0.0868   0.0458

sigma^2 = 6300:  log likelihood = -5907.95
AIC=11833.91   AICc=11834.09   BIC=11878.26

Training set error measures:
                  ME     RMSE      MAE MPE MAPE      MASE        ACF1
Training set 0.44309 79.02328 45.78169 NaN  Inf 0.5805454 0.003263827

Code

## Forecasting Test Data -------------------------------------------------------
east_values <- forecast(east_arima, h = length(east_test_ts))
midwest_values <- forecast(midwest_arima, h = length(midwest_test_ts))
south_values <- forecast(south_arima, h = length(south_test_ts))
west_values <- forecast(west_arima, h = length(west_test_ts))

8.2 Accuracy

Code

## Plotting & Accuracy

## East Region
east_auto_acc <- forecast::accuracy(east_values)

east_auto_acc  %>% 
  kbl(caption = "<center>Auto ARIMA East Accuracy<center>") %>% 
  kable_styling()

Auto ARIMA East Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.1506172	82.50546	47.1043	NaN	Inf	0.7224607	0.0170736

Code

east_auto_plot <- forecast::autoplot(east_values) +
  geom_line(data = east_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "East Region Auto ARIMA",
       y = "Daily Deaths")

east_auto_plot 

Code

## Midwest Region
midwest_auto_acc <- forecast::accuracy(midwest_values)

midwest_auto_acc %>% 
  kbl(caption = "<center>Auto ARIMA Midwest Accuracy<center>") %>% 
  kable_styling()

Auto ARIMA Midwest Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.385353	102.0896	59.56282	NaN	Inf	0.5635082	-0.0553299

Code

midwest_auto_plot <- forecast::autoplot(midwest_values) +
  geom_line(data = midwest_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "Midwest Region Auto ARIMA",
       y = "Daily Deaths")

midwest_auto_plot 

Code

## South Region
south_auto_acc <- forecast::accuracy(south_values)

south_auto_acc %>% 
  kbl(caption = "<center>Auto ARIMA South Accuracy<center>") %>% 
  kable_styling()

Auto ARIMA South Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.7155848	206.0793	126.7407	NaN	Inf	0.7335323	0.0086942

Code

south_auto_plot <- forecast::autoplot(south_values) +
  geom_line(data = south_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "South Region Auto ARIMA",
       y = "Daily Deaths")

south_auto_plot

Code

## West Region
west_auto_acc <- forecast::accuracy(west_values)

west_auto_acc %>% 
  kbl(caption = "<center>Auto ARIMA West Accuracy<center>") %>% 
  kable_styling()

Auto ARIMA West Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.44309	79.02328	45.78169	NaN	Inf	0.5805454	0.0032638

Code

west_auto_plot <- forecast::autoplot(west_values) +
  geom_line(data = west_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "West Region Auto ARIMA",
       y = "Daily Deaths")

west_auto_plot

Code

# saving results
save(east_auto_plot,
     midwest_auto_plot,
     south_auto_plot,
     west_auto_plot,
     east_auto_acc,
     midwest_auto_acc,
     south_auto_acc,
     west_auto_acc,
     file = "results/uni_autoarima.rda")

9. Compare Arima Models

Code

load("results/uni_arima.rda")
load("results/uni_sarima.rda")
load("results/uni_autoarima.rda")

East Results

Code

east1 <- east_initial_arima_acc %>% 
  as.data.frame()

east2 <- east_tuned_arima_acc %>% 
  as.data.frame()

east3 <- east_initial_sarima_acc%>% 
  as.data.frame() 

east4 <- east_tuned_sarima_acc%>% 
  as.data.frame()

east5 <- east_auto_acc%>% 
  as.data.frame()

tribble(
  ~model, ~mae, ~mase,
  #--|--|----
  "initial arima", east1$MAE, east1$MASE,
  "tuned arima", east2$MAE, east2$MASE,
  "initial sarima", east3$MAE, east3$MASE,
  "tuned sarima", east4$MAE, east4$MASE,
  "auto arima", east5$MAE, east5$MASE,
) %>% 
  kbl() %>% 
  kable_styling()

model	mae	mase
initial arima	64.66739	0.9918341
tuned arima	47.18823	0.7237480
initial sarima	61.27670	0.9398295
tuned sarima	53.94287	0.8273471
auto arima	47.10430	0.7224607

Midwest Results

Code

midwest1 <- midwest_initial_arima_acc %>% 
  as.data.frame()

midwest2 <- midwest_tuned_arima_acc %>% 
  as.data.frame()

midwest3 <- midwest_initial_sarima_acc%>% 
  as.data.frame() 

midwest4 <- midwest_tuned_sarima_acc%>% 
  as.data.frame()

midwest5 <- midwest_auto_acc%>% 
  as.data.frame()

tribble(
  ~model, ~mae, ~mase,
  #--|--|----
  "initial arima", midwest1$MAE, midwest1$MASE,
  "tuned arima", midwest2$MAE, midwest2$MASE,
  "initial sarima", midwest3$MAE, midwest3$MASE,
  "tuned sarima", midwest4$MAE, midwest4$MASE,
  "auto arima", midwest5$MAE, midwest5$MASE,
) %>% 
  kbl() %>% 
  kable_styling()

model	mae	mase
initial arima	78.32550	0.7410171
tuned arima	54.43655	0.5150099
initial sarima	78.32550	0.7410171
tuned sarima	63.95655	0.6050762
auto arima	59.56282	0.5635082

South Results

Code

south1 <- south_initial_arima_acc %>% 
  as.data.frame()

south2 <- south_tuned_arima_acc %>% 
  as.data.frame()

south3 <- south_initial_sarima_acc%>% 
  as.data.frame() 

south4 <- south_tuned_sarima_acc%>% 
  as.data.frame()

south5 <- south_auto_acc%>% 
  as.data.frame()

tribble(
  ~model, ~mae, ~mase,
  #--|--|----
  "initial arima", south1$MAE, south1$MASE,
  "tuned arima", south2$MAE, south2$MASE,
  "initial sarima", south3$MAE, south3$MASE,
  "tuned sarima", south4$MAE, south4$MASE,
  "auto arima", south5$MAE, south5$MASE,
) %>% 
  kbl() %>% 
  kable_styling()

model	mae	mase
initial arima	129.55599	0.7498262
tuned arima	94.70831	0.5481396
initial sarima	129.55599	0.7498262
tuned sarima	94.70831	0.5481396
auto arima	126.74072	0.7335323

West Results

Code

west1 <- west_initial_arima_acc %>% 
  as.data.frame()

west2 <- west_tuned_arima_acc %>% 
  as.data.frame()

west3 <- west_initial_sarima_acc%>% 
  as.data.frame() 

west4 <- west_tuned_sarima_acc%>% 
  as.data.frame()

west5 <- west_auto_acc%>% 
  as.data.frame()

tribble(
  ~model, ~mae, ~mase,
  #--|--|----
  "initial arima", west1$MAE, west1$MASE,
  "tuned arima", west2$MAE, west2$MASE,
  "initial sarima", west3$MAE, west3$MASE,
  "tuned sarima", west4$MAE, west4$MASE,
  "auto arima", west5$MAE, west5$MASE,
) %>% 
  kbl() %>% 
  kable_styling()

model	mae	mase
initial arima	62.19197	0.7886397
tuned arima	46.88033	0.5944769
initial sarima	62.19197	0.7886397
tuned sarima	55.78604	0.7074078
auto arima	45.78169	0.5805454